﻿/////////////////////////////////////////////////////////////////////
//
//	Map_FVG C# Class Library.
//
//	UTM
//	conversione coordinate geodetiche in utm e viceversa
//  Documenti/Mappe Geografiche/Geodesics/GeoTrans/win_dev/geotrans3.3/CCS/sc/dtcc/CoordinateSystems/transmerc
//
//	autore: Gustavo Mainardi 
//	versione: 1.0
//	data: 01/12/2015
//
//
//	Cronologia delle versioni:
//
//	Versione 1.0 01/01/2015
//		versione originale
//
//////////////////////////////////////////////////////////////////////////
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
namespace Map_Fvg
{
    public partial class UTM 
    {
        double PI, MAX_LAT, MAX_DELTA_LONG, PI_OVER_180;

        double TranMerc_es;             /* Eccentricity squared */
        double TranMerc_ebs;            /* Second Eccentricity squared */

        /* Transverse_Mercator projection Parameters */
        double TranMerc_Origin_Lat;           /* Latitude of origin in radians */
        double TranMerc_Origin_Long;          /* Longitude of origin in radians */
        double TranMerc_False_Northing;       /* False northing in meters */
        double TranMerc_False_Easting;        /* False easting in meters */
        double TranMerc_Scale_Factor;         /* Scale factor  */

        /* Isometric to geodetic latitude parameters */
        double TranMerc_ap;
        double TranMerc_bp;
        double TranMerc_cp;
        double TranMerc_dp;
        double TranMerc_ep;

        double semiMajorAxis;

        public struct UTMpoint
        {
            public double lat;
            public double lon;
            public int easting;
            public int northing;
            public int gb_X;
            public int gb_Y;
        }

        UTMpoint geodeticCoordinates;

         public UTM()
        {
   //         InitializeComponent();

         }

        public void geoToUTM()
        {
            PI      =     3.14159265358979323e0;
            MAX_LAT    =  (( 84.5 * PI) / 180.0); /*  84.5 degrees in radians */
            PI_OVER_180 = (3.14159265358979323e0 / 180.0);
            double PI_OVER_2    =     (PI/2.0);
            MAX_DELTA_LONG  =  ((PI * 90)/180.0);
            double ellipsoidSemiMajorAxis = 6378137.0;
            double ellipsoidFlattening =  1 / 298.257223563;
            double inv_f = 1 / ellipsoidFlattening;

            double centralMeridian = 15.0*PI_OVER_180;      //  meridiano centrale per il FVG
            double originLatitude = 0;  
            double falseEasting   = 500000;
            double falseNorthing  = 0;
            double scaleFactor = 0.9996;

            semiMajorAxis = ellipsoidSemiMajorAxis;
            double flattening = ellipsoidFlattening;

            double tn;         /* True Meridianal distance constant  */
            double tn2;
            double tn3;
            double tn4;
            double tn5;

            TranMerc_Origin_Lat = originLatitude;
            TranMerc_Origin_Long    = centralMeridian;
            TranMerc_False_Northing = falseNorthing;
            TranMerc_False_Easting  = falseEasting; 
            TranMerc_Scale_Factor   = scaleFactor;

  /* Eccentrie city Squared */
            TranMerc_es = 2 * flattening - flattening * flattening;

  /* Second Eccentricity Squared */
            TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;

/* Semi-minor axis of ellipsoid, in meters */
            double TranMerc_b = semiMajorAxis * (1 - flattening); 
   
  /*True meridianal constants  */
            tn = (semiMajorAxis - TranMerc_b) / (semiMajorAxis + TranMerc_b);
            tn2 = tn * tn;
            tn3 = tn2 * tn;
            tn4 = tn3 * tn;
            tn5 = tn4 * tn;

            TranMerc_ap = semiMajorAxis * 
                (1.0 - tn + 5.0 * (tn2 - tn3)/4.0 + 81.0 * (tn4 - tn5)/64.0 );
            TranMerc_bp = 3.0 * semiMajorAxis * 
               (tn - tn2 + 7.0 * (tn3 - tn4) /8.0 + 55.0 * tn5/64.0 ) / 2.0;
            TranMerc_cp = 15.0 * semiMajorAxis *
                (tn2 - tn3 + 3.0 * (tn4 - tn5 )/4.0) /16.0;
            TranMerc_dp = 35.0 * semiMajorAxis *
                (tn3 - tn4 + 11.0 * tn5 / 16.0) / 48.0;
            TranMerc_ep = 315.0 * semiMajorAxis * (tn4 - tn5) / 512.0;

            convertToWGS84();
            convertToGaussBoaga();

        }

  
        int convertFromGeodetic( )
        {
/*
 * The function convertFromGeodetic converts geodetic
 * (latitude and longitude) coordinates to Transverse Mercator projection
 * (easting and northing) coordinates, according to the current ellipsoid
 * and Transverse Mercator projection coordinates.  If any errors occur, 
 * an exception is thrown with a description of the error.
 *
 *    longitude     : Longitude in radians                        (input)
 *    latitude      : Latitude in radians                         (input)
 *    easting       : Easting/X in meters                         (output)
 *    northing      : Northing/Y in meters                        (output)
 */

            double c;       /* Cosine of latitude                          */
            double c2;
            double c3;
            double c5;
            double c7;
            double dlam;    /* Delta longitude - Difference in Longitude       */
            double eta;     /* constant - TranMerc_ebs *c *c                   */
            double eta2;
            double eta3;
            double eta4;
            double s;       /* Sine of latitude                        */
            double sn;      /* Radius of curvature in the prime vertical       */
            double t;       /* Tangent of latitude                             */
            double tan2;
            double tan3;
            double tan4;
            double tan5;
            double tan6;
            double t1;      /* Term in coordinate conversion formula - GP to Y */
            double t2;      /* Term in coordinate conversion formula - GP to Y */
            double t3;      /* Term in coordinate conversion formula - GP to Y */
            double t4;      /* Term in coordinate conversion formula - GP to Y */
            double t5;      /* Term in coordinate conversion formula - GP to Y */
            double t6;      /* Term in coordinate conversion formula - GP to Y */
            double t7;      /* Term in coordinate conversion formula - GP to Y */
            double t8;      /* Term in coordinate conversion formula - GP to Y */
            double t9;      /* Term in coordinate conversion formula - GP to Y */
            double tmd;     /* True Meridianal distance                        */
            double tmdo;    /* True Meridianal distance for latitude of origin */
            double temp_Origin;
            double temp_Long;

            double longitude = geodeticCoordinates.lon * PI_OVER_180;
            double latitude = geodeticCoordinates.lat * PI_OVER_180;
           
            if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
            {  /* Latitude out of range */
                MessageBox.Show("valore anomalo di latitudinte:  " + latitude.ToString() );
                return 1;
            }

            if (longitude > PI) longitude -= (2 * PI);
            if ((longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
                 || (longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
            {
                 if (longitude < 0)
                     temp_Long = longitude + 2 * PI;
                else
                    temp_Long = longitude;
                if (TranMerc_Origin_Long < 0)
                    temp_Origin = TranMerc_Origin_Long + 2 * PI;
                else
                    temp_Origin = TranMerc_Origin_Long;
                if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))|| (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
                {
                    longitude = geodeticCoordinates.lon;
                    MessageBox.Show("valore anomalo di longitudine:  " + latitude.ToString() );
                    return 1;
                }
            }

  /* 
   *  Delta Longitude
   */
            dlam = longitude - TranMerc_Origin_Long;

            if (Math.Abs(dlam) > (9.0 * PI_OVER_180))
            { /* Distortion results if Longitude is > 9 delongituderees from Central Meridian */
                MessageBox.Show("warning: longitude " + longitude.ToString() );
    //            return 1;  //  ?
            }

            if (dlam > PI) dlam -= (2 * PI);
            if (dlam < -PI) dlam += (2 * PI);
            if (Math.Abs(dlam) < 2.0-10) dlam = 0.0;

            s  = Math.Sin(latitude);
            c = Math.Cos(latitude);
            c2 = c * c;
            c3 = c2 * c;
            c5 = c3 * c2;
            c7 = c5 * c2;
            t = Math.Tan(latitude);
            tan2 = t * t;
            tan3 = tan2 * t;
            tan4 = tan3 * t;
            tan5 = tan4 * t;
            tan6 = tan5 * t;
            eta  = TranMerc_ebs * c2;
            eta2 = eta * eta;
            eta3 = eta2 * eta;
            eta4 = eta3 * eta;

  /* radius of curvature in prime vertical */
            sn  = sphsn(s);

  /* True Meridianal Distances */
            tmd = sphtmd(latitude, s, c);

  /*  Origin  */
            tmdo = sphtmd(TranMerc_Origin_Lat, Math.Sin(TranMerc_Origin_Lat), Math.Cos(TranMerc_Origin_Lat));

  /* northing */
            t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
            t2 = sn * s * c * TranMerc_Scale_Factor/ 2.0;
            t3 = sn * s * c3 * TranMerc_Scale_Factor * 
                 (5.0 - tan2 + 9.0 * eta + 4.0 * eta2) /24.0; 

             t4 = sn * s * c5 * TranMerc_Scale_Factor * (61.0 - 58.0 * tan2
                + tan4 + 270.0 * eta - 330.0 * tan2 * eta + 445.0 * eta2
                + 324.0 * eta3 -680.0 * tan2 * eta2 + 88.0 * eta4 
                -600.0 * tan2 * eta3 - 192.0 * tan2 * eta4) / 720.0;

            t5 = sn * s * c7 * TranMerc_Scale_Factor *
                (1385.0 - 3111.0 * tan2 + 543.0 * tan4 - tan6) / 40320.0;

            double dlamSq = dlam * dlam;

            double northing = TranMerc_False_Northing + t1 + 
                (((t5 * dlamSq + t4) * dlamSq + t3) * dlamSq + t2) * dlamSq;

  /* Easting */
            t6 = sn * c  * TranMerc_Scale_Factor;
            t7 = sn * c3 * TranMerc_Scale_Factor * (1.0 - tan2 + eta ) /6.0;
            t8 = sn * c5 * TranMerc_Scale_Factor * (
                5.0 - 18.0 * tan2 + tan4 + 14.0 * eta - 58.0 * tan2 * eta +
                13.0 * eta2 + 4.0 * eta3 - 64.0 * tan2 * eta2 - 24.0 * tan2 * eta3 ) / 120.0;
            t9 = sn * c7 * TranMerc_Scale_Factor *
                ( 61.0 - 479.0 * tan2 + 179.0 * tan4 - tan6 ) /5040.0;

            double easting = TranMerc_False_Easting + 
                (((t9 * dlamSq + t8) * dlamSq + t7) * dlamSq + t6) * dlam;

            geodeticCoordinates.easting = (int)easting;
            geodeticCoordinates.northing = (int)northing;

           return 0;

        }

        double sphsn( double  sinLat ) 
        {
        //   double sinLat = Math.sin(latitude);
            return semiMajorAxis / Math.Sqrt( 1.0 - TranMerc_es * sinLat * sinLat);
        }


        double sphsr( double sinLat) 
        {
            double denom = Math.Sqrt(1.0 - TranMerc_es * sinLat * sinLat);
            return semiMajorAxis * (1.0 - TranMerc_es) / denom / denom / denom;
        }

        double sphtmd( double latitude, double sinLat, double cosLat )
        {
            double sin2Lat = 2.0 * sinLat * cosLat;
            double cos2Lat = 2.0 * cosLat * cosLat - 1.0;
            double sin2LatCubed = sin2Lat * sin2Lat * sin2Lat;
            double sin4Lat = 2.0 * sin2Lat * cos2Lat;
            double sin6Lat = 3.0 * sin2Lat - 4.0 * sin2LatCubed;
            double sin8Lat = cos2Lat * ( 4.0 * sin2Lat - 8.0 * sin2LatCubed);

            return TranMerc_ap * latitude - TranMerc_bp * sin2Lat + TranMerc_cp * sin4Lat
                        - TranMerc_dp * sin6Lat + TranMerc_ep * sin8Lat;

        }

        void convertToWGS84()
        {
            for (int i = 0; i < GPS.nPoints; ++i)
            {
               geodeticCoordinates.lon = GPS.vectG[i].lon;
               geodeticCoordinates.lat = GPS.vectG[i].lat;

                int iret = convertFromGeodetic();
                if (iret > 0)
                {

                    MessageBox.Show("errore: convertFromGeodetic " + iret.ToString());
                    return;
                }
                GPS.vectG[i].easting = geodeticCoordinates.easting;
                GPS.vectG[i].northing = geodeticCoordinates.northing;

            }
        }

        void convertToGaussBoaga()
        {
            int shift_X = 2020004;
            int shift_Y = 23;
            for (int i = 0; i < GPS.nPoints; ++i)
            {
                GPS.vectG[i].gb_X = GPS.vectG[i].easting + shift_X;
                GPS.vectG[i].gb_Y = GPS.vectG[i].northing + shift_Y;
            }
        }

        int convertFromUTM()
        {

            return (0);
        }

    }
}
